# Copyright (c) HySoP 2011-2024
#
# This file is part of HySoP software.
# See "https://particle_methods.gricad-pages.univ-grenoble-alpes.fr/hysop-doc/"
# for further info.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import sympy as sm
from hysop.constants import __VERBOSE__, __DEBUG__
from hysop.tools.henum import EnumFactory
from hysop.tools.htypes import check_instance
from hysop.numerics.remesh.kernel_generator import Kernel, SymmetricKernelGenerator
Remesh = EnumFactory.create(
"Remesh",
[
"L1_0",
"L2_1",
"L2_2",
"L4_2",
"L4_4",
"L6_4",
"L6_6",
"L8_4", # lambda remesh kernels
"L2_1s",
"L2_2s",
"L4_2s",
"L4_4s",
"L6_4s",
"L6_6s",
"L8_4s", # splitted lambda remesh kernels
"Mp4",
"Mp6",
"Mp8", # Mprimes kernels: Mp4 = M'4 = L2_1 and Mp6 = M'6 = L4_2
"M4",
"M8", # M kernels
"O2",
"O4", # Corrected kernels, allow a large CFL number
"L2", # Corrected and limited lambda 2
],
)
[docs]
class RemeshKernelGenerator(SymmetricKernelGenerator):
[docs]
class RemeshKernel(Kernel):
def __init__(
self,
moments,
regularity,
verbose=__DEBUG__,
split_polys=False,
override_cache=False,
):
generator = RemeshKernelGenerator(verbose=verbose)
generator.configure(n=moments)
kargs = generator.solve(
r=regularity,
override_cache=override_cache,
split_polys=split_polys,
no_wrap=True,
)
super().__init__(**kargs)
[docs]
@staticmethod
def from_enum(remesh):
check_instance(remesh, Remesh)
remesh = str(remesh)
assert (
remesh[0] == "L" and (remesh != "L2") or (remesh in ("M4", "M8"))
), "Only lambda remesh kernels are supported."
if remesh in ("M4", "M8"):
# given M4 or M8 kernels
x = sm.Symbol("x")
if remesh == "M4":
M4 = (
sm.Poly(
(1 / sm.Rational(6)) * ((2 - x) ** 3 - 4 * (1 - x) ** 3), x
),
sm.Poly((1 / sm.Rational(6)) * ((2 - x) ** 3), x),
)
return Kernel(
n=2,
r=4,
deg=3,
Ms=2,
Mh=None,
H=None,
remesh=True,
P=(M4[1].subs(x, -x), M4[0].subs(x, -x), M4[0], M4[1]),
)
else:
remesh = remesh[1:]
if remesh[-1] == "s":
remesh = remesh[:-1]
split_polys = True
else:
split_polys = False
remesh = [int(x) for x in remesh.split("_")]
assert len(remesh) == 2
assert remesh[0] >= 1
assert remesh[1] >= 0
return RemeshKernel(remesh[0], remesh[1], split_polys=split_polys)
def __str__(self):
return f"RemeshKernel(n={self.n}, r={self.r}, split={self.poly_splitted})"
if __name__ == "__main__":
import numpy as np
from matplotlib import pyplot as plt
for i in range(1, 5):
p = 2 * i
kernels = []
for r in [1, 2, 4, 8]:
try:
kernel = RemeshKernel(p, r)
kernels.append(kernel)
except RuntimeError:
print(f"Solver failed for p={p} and r={r}.")
if len(kernels) == 0:
continue
k0 = kernels[0]
fig = plt.figure()
plt.xlabel(r"$x$")
plt.ylabel(r"$\Lambda_{" + "{},{}".format(p, "r") + "}$")
X = np.linspace(-k0.Ms - 1, +k0.Ms + 1, 1000)
s = plt.subplot(1, 1, 1)
for i, k in enumerate(kernels):
s.plot(X, k(X), label=r"$\Lambda_{" + f"{p},{k.r}" + "}$")
s.plot(k0.I, k0.H, "or")
axe_scaling = 0.10
ylim = s.get_ylim()
Ly = ylim[1] - ylim[0]
s.set_ylim(ylim[0] - axe_scaling * Ly, ylim[1] + axe_scaling * Ly)
s.legend()
plt.show(block=True)
# plt.savefig('out/gamma_{}_r'.format(p))